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i.  mrsowcnat  and  summary. 

This  paper  describes  a  FORTRAN  subroutine  called  WNEW,  for  approximating  any 
one  of  the  following  four  Integrals: 


(1.1) 

f  f  (x)dx  ; 

*  *00 

(1.2) 

1  f(t)dt  (f  not  oscillatory,  A  finite) 
>  A 

(1.3) 

j  F(t)dt  (F  oscillatory,  A  finite) 

'A 

rB 

(1.4) 

f(t)dt  (A,B  finite) 

J  A 

We  remark  in  view  of  (1.2)  and  (1.3)  that 


rB 


f(t)dt 


A  description  of  the  parameters  of  the  subroutine  WNEW  and  the  method  of 
calling  it  are  given  in  Sec.  6  of  this  paper;  the  user  vho  does  not  wish  to 
concern  himself  with  the  special  powers  or  pitfalls  of  this  subroutine  should 
skip  directly  to  Sec^6. 

The  formulas  that  this  subroutine  is  based  on  are  most  powerful  when  the 

A 

integrand  f  does  not  have  a  singularity  (a  singularity  is  a  point  there 
df/dx  does  not  exist)  in  the  interior  of  the  range  of  integration;  however, 
singularities  at  end-points  of  intervals  are  allowed.  Indeed,  it  is  in  the 
cases  where  f  has  singularities  at  the  end-points  of  the  range  of  integration^/ 
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that  the  subroutine  WNEW  is  superior  to  other  subroutines . concede  that 
Gaussian  quadrature  which  houses  the  singularities  of  the  integrand  in  the 
weight  function  [1,2]  may  be  superior  to  the  WNEW  methods .  However ,  no  method  [8] 
is  superior  to  the  WNEW  method  if  the  exact  nature  of  the  singularities  at 
end  points  of  intervals  are  unknown,  or  ignored.  The  subroutines  may  then 
be  used  either  directly,  or  be  a  part  of  a  polyalgorithn  [4,5]  which  gets 
accurate  results  in  spite  of  singularities  of  the  integrand  at  end-points 
of  an  interval  of  integration. 

We  emphasize  that  the  WNEW  subroutine  yields  accurate  results  in  spite 
of  singularities  at  end  points  of  an  interval.  The  presence  of  singularities 
at  interior  points  of  the  interval  of  integrations  any  considerably  slow  up 
the  rate  of  convergence.  If  the  function  f  in  the  integral 

rd 

(1.5)  f (t)dt 

Jc 

(where  c  or  d  may  be  either  finite  or  infinite)  has  singularities  in  the 
interval  (c,d)  ,  then  in  order  to  achieve  best  accuracy  we  strongly  recommend 
replacing  (1.5)  by  a  finite  number  of  integrals  with  the  property  that  each 
only  has  singularities  at  the  end-points  of  an  interval.  For  example,  if  the 
function  f  in  (1.5)  has  singularities  at  u  and  v  ,  vrtiere  c  <  u  <  v  <  d  , 
then  we  reconmend  replacing  (1.5)  as  follows: 

/d  rU  fV  rd 

(1.6)  f(t)dt  =  f(t)dt  +  f(t)dt  +  f(t)dt  . 

'  C  'C  ■'u  'V 

Each  of  these  integrals  may  now  be  accurately  approximated  by  an  appropriate 
formula  used  to  approximate  (1.2),  (1.3)  or  (1.4). 


The  subroutine  WNEW  is  based  on  the  theory  of  [6,7  J ;  see  also  the  summary 
paper  [3].  In  Sec.  2  we  briefly  surmarize  the  transformation  used  to  transform 
each  of  the  integrals  (1.2),  (1.3)  or  (1.4)  into  (1.1),  the  error  when  the 
trapezoidal  formula  is  applied  to  (1.1),  as  well  as  a  more  accurate  description 
of  the  type  of  integrals  for  which  the  formulas  for  approximating  (1.1),  (1.2), 
(1.3),  and  (1.4)  are  most  effective. 

Section  3  describes  the  basis  of  the  algorithm  by  conbining  trapezoidal 
and  midordinate  rules .  In  Sec .  4  we  give  some  examples  which  illustrate  the 
application  of  the  algorithm.  In  Sec.  5  we  illustrate  some  pitfalls  of  the 
algorithm,  arising  as  a  consequence  of  inaccurately  computing  the  integrand 
near  a  singularity.  We  also  illustrate  methods  of  circumventing  these  pitfalls. 

In  Sec.  6  we  give  a  precise  description  of  the  subroutine 

WOT(INTRUL,A,B,EPS,IP) 

and  of  the  role  of  the  parameters  in  this  subroutine.  We  also  give  a  flow¬ 
chart  description  of  the  main  ideas  of  the  subroutine. 

In  Sec.  7  we  give  an  explicit  FORTRAN  listing  of  the  subroutine  VCEW. 
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2.  BASIC  IDEAS,  TRANSFORMATIONS  AND  ERROR. 

The  algorithms  of  the  program  are  all  based  on  the  trapezoidal  formula 

fOO  CO 

(2.1)  [  f  (x)dx  ^  h  l  f(kh) 

J  —CO  lc=-oo 

where  h  >  0  is  the  step  size.  If  f  has  no  singularities  cn  (-«>,<»)  (e.g. 

if  a  single  formula  is  used  to  describe  f  on  the  whole  interval  (-00  ,00)) 
then  it  may  be  shown  that  the  error  of  formula  (2.1)  satisfies 

(2.2)  | error.  <  C  e"c/h 

where  C  and  c  are  positive  constants  that  are  independent  of  h  .  Thus  if 
h  is  replaced  by  h/2  ,  then  the  correct  number  of  significant  figures  in 
the  approximation  (2.1)  doubles.  Best  results  are  achieved  for  (2.1)  if  in 
addition  to  being  analytic  on  1  ,  f  also  satisfies  the  inequality 

(2.3)  |f(x)j  <  C,e“aM 

on  the  real  line  TR  ,  where  C'  and  a  are  positive  constants.  In  this  case 
relatively  few  points  are  required  in  the  trapezoidal  sun  to  achieve  the 
desired  accuracy.  If  f  decreases  to  zero  at  an  algebraic  rate,  as  t  -*•  ±«®  , 
the  formula  (2.1)  is  still  accurate,  however  in  that  case  many  more  points 
are  required  to  achieve  a  desired  accuracy.  This  latter  situation  can  sometimes 
be  remedied  by  use  of  the  transformation 


I 


•  ^  V 

► 


(2.4) 


t  =  xe 


and  then  applying  the  trapezoidal  formula  to  the  transformed  integral. 

The  integral  (1.2)  is  transformed  into  the  integral  (1.1)  by  means  of  the 
transformation 


(2.5) 


t  «  A  +  e 


Applying  (2.1)  after  making  the  transformation  (2.5)  results  in  the  quadrature 


formula 


(2.6) 


f(t)dt  sh  i  ekhf(A+ekh) 

■  A 


(We  remark  here,  that 


(2.7)  [B  f(t)dt  =  f  f(-t)dt  s  h  [  e^hf(B-e^th)  .) 

J  -oo  J  -B  k=-°° 

The  formula  (2.6)  is  particularly  accurate  when  used  to  approximate  integrals 
for  which  the  integrands  have  an  algebraic-type  singularity  at  A  ,  and  vrtiich 
approach  zero  at  an  algebraic  rate  as  x  -*■  »  ,  such  as  integrals  of  the  form 

°°  t’^t  +  D^dt  ,  or  [  t3"r(logt)(l+t2)"%dt  . 

^0  J0 

If  A  =  0  ,  the  ideal  boundedness  condition  on  f  corresponding  to  (2.3)  is 


C’t01"1  ,  0  <  t  <  1 


!f(t)|<  _!_a 

ICt  a  ,  t  >  1  ; 


(2.8) 
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if  (2.3)  is  satisfied  then  after  making  the  transformation  (2.5)  (with  A  =  0) 
one  gets  an  integral  over  ]R  for  which  the  integrand  satisfies  (2.3). 

The  integral  (1.3)  is  transformed  into  the  integral  (1.1)  by  means  of  the 
transformation 

(2.9)  t  =  A  +  log{ex  +  (l+e2x)%)  . 

After  making  the  transformation  (2.9)  and  applying  (2.1),  we  get* 

(2.10)  '  f(t)dt  =  h£  (l  +  e"2kh)‘%f(A+  lpgte1*  + /l  +  e21*  })  . 

JA  k=-“ 

This  formula  is  best  suited  for  the  evaluation  of  integrals  for  which  the 
integrand  f(t)  has  an  algebraic-type  singularity  at  t  -  A  and  which  be¬ 
haves  in  an  oscillatory  manner  as  t  •*  ®  .  Examples  of  such  integrals  are 

foo  .  .  .  ,  .  . 

t”^7  c”cos(3t)dt  ,  or  log[l  -  — ']e~  at  .  If  A  -  0  in  (1.3), 

Jo  Jo 

the  ideal  boundedness  condition  on  f  corresponding  to  (2.3)  is 
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if  f  satisfies  (2.11)  then  after  making  the  transformation  (2.9)  (with  A  =  0) 
one  gets  an  integral  over  H  for  which  the  integrand  satisfies  (2.3).  If  f 
is  oscillatory  on  (0,°°)  but  does  not  decrease  to  zero  at  the  rate  (2.11),  as 
in  the  case  of  the  evaluation  of  some  semi- infinite  transforms,  such  as 

|°°  t"^jQ(at)dt  the  formula  (2.10)  is  still  quite  accurate.  However,  in  this 

case  a  large  number  of  points  are  required  to  achieve  a  desired  accuracy.  This 
situation  may  be  remedied  *  by  an  Euler  technique,  such  as  that  described  in  [1]. 

The  integral  (1.4)  is  transformed  into  the  integral  (1.1)  by  means  of  the 
transforms  ion 


(2.12) 


t  = 


A+Bex 
l  +  ex 


The  boundedness  condition  on  f  corresponding  to  (2.3)  is 


(2.13)  |f(t)|  <  C|(t-A)(B-t)r"1  ,  A  <  t  <  B 


where  C  and  a  are  positive  constants.  Making  the  transformation  (2.12) 
and  then  applying  the  trapezoidal  formula,  we  get 

fB  r  e1*  A+Be^1 

(2.14)  f(t)dt  *  (B-A)hJ  f^TTRK)  • 

;A  k=-»  (1  +  e  )  1  +  e 


^Another  method  of  circumventing  this  difficulty  is  to  evaluate 

,00 

10)  =  |^e"Xtt-1^I0(at)dt  for  e.g.  A  =  1/2, 1/4  and  1/8  ,  and  then  extrapolate 
to  the  limit  A  =  0  . 
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3.  BASIS  OF  ALGORITHM. 

Let  us  define  T^(f)  and  ll^(f)  by 

T.  (f)  •hi  f  (kh) 

h  k=-°° 

(3.1) 

M(f)  =  h  l  f((2k-l)h/2)  . 

k=-oo 

Thus  the  sum  on  the  right  hand  side  of  (2.1)  is  T^(f)  .  Furthermore,  it 
follows  that 

(3.2)  Th/2(f)  =  %[Th(f)  +Mh(f)1  • 

Let  us  start  with  h  =  1  (say)  and  then  compute  T^(f)  .  The  bound  (2.2) 
shows  that  has  at  least  twice  as  many  significant  figures  as  T^(f) 

Next,  let  us  compute  M^(f)  ,  as  well  as  the  difference 

(3.3)  |Th(f)-K  (f)  |  =  c/3  Ul\\  f(x)dx-Th(f)]l  . 

U  -00  ^  ' 

Thus,  for  sufficiently  small  e  , 

(3.4)  |f  f(x)dx-Th/2(f)|  =  |f  f(x)dx-%[Th(f)+Nh(f)]|  =  0(e2)  <  c  . 

JC\  ^Ju 

In  practice,  we  cannot  sum  all  of  the  terms  in  the  infinite  suns  (3.1). 
The  assumption  (2.3)  then  offers  a  convenient  stopping  criteria.  Suppose  that 
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we  stop  the  s urination  in  T^(f)  for  k  >  0  when 
(3.5)  |f(Nh)|(  =  0(e‘aNh))  <  e/3  . 

Then  we  may  expect  that 


(3.6) 


h  l  f(kh) 
k-Nfl 


<  ofh  i  e'1*] 


=  0 


rh  g-aCNhDh'i 


t  l-e"ah  > 


=  0(e_aNh)  =  0(e) 


That  is,  we  may  expect  the  tail  of  the  series  to  be  of  the  same  order  of 
magnitude  as  last  included  term.  In  order  to  avoid  stopping  the  algorithm 
at  or  near  a  zero  of  f  in  practice,  we  make  the  more  reliable  test 

(3.7)  |f  (Nh)  |  +  |f((Nfl)h)|  +  Jf(.(W-2)h)|  <  |  . 

Similarly,  the  sums  on  the  right  hand  sides  of  (2.6),  (2.10)  and  (2.14) 
share  the  properties  of  the  trapezoidal  formula  on  the  right  hand  side  of  (2.1), 
under  the  assumptions  of  (2.8),  (2.11)  and  (2.13)  respectively,  which  correspond 
to  the  assumption  (2.3). 
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4.  EXILES. 

The  examples  of  this  section  illustrate  the  applications  of  each  of  the 
formulas  (2.1),  (2.6),  (2.10)  and  (2.14).  Each  of  these  formulas  may  be 
represented  as  a  single  approximating  formula 

f  L 

(4.1)  I  f(x)dx  s  h  £  w^(h)f (^(h)) 

i  k=-K 

The  results  of  various  examples  are  tabulated  in  Table  4.1.  In  this  table  we 
tabulate  the  integral  to  be  approximated,  the  exact  value  of  the  integral,  the 
method  of  quadrature  used,  EPS,  a  parameter  specified  by  the  user  and  which 
is  the  s  of  the  previous  section,  K,L  (see  Eq.  (4.1))  and  the  final 
approximation  achieved. 

In  the  third  last  and  second  last  entry  of  Table  4.1  it  was  not  possible 
to  achieve  the  accuracy  of  e  .  This  phenomenon  occurs  due  to  a  pitfall  in 
computations ;  such  pitfalls  and  ways  of  circumventing  these  are  explained  in 
the  next  section. 


t 


TABLE  4.1  EXAMPLES 


{  INTEGRAL 

METHOD 

USED* 

EXACT 

ANSWER 

EPS 

K 

B 

ABSOLUTE  VALUE" 
ERROR  M 

FINAL  APPKOXIMAX] 

*• 

$ 

r00  2 

e  x  dx 

— oo 

Eq.  (2.1) 

/if 

10'7 

10'16 

18 

36 

18 

36 

<  10'7 
<  10-15 

2 

4r- 

o  4* 

J  ** 

Eq.  (2.14) 

TT  = 

.  10'7 
10'9 

78 

180 

42 

108 

<  10-7 
<  10-16 

f  x  Vxdx 

0 

Eq.  (2.10) 

/f 

10-16 

648 

312 

<  10'16 

r°°  v  % 

l0  1  +  X  * 

Eq.  (2.6) 

*r 

10-11 

114 

114 

<  10'12 

J 

log[l-s“(lVXdx 

0  x 

Eq.  (2.10) 

-3.045689266 

50352 

10~14 

168 

132 

<  10-16 

hi 

f 

(3-2x-xZ)  ^dx 

-1 

Eq.  (2.14) 

7T 

10~3 

10"7 

12 

21 

18 

36 

- - - -j' 

_A  i 

<  Hf*i 

| 

.  2 

•2  dx 

0  [x(4-x)]^ 

Eq.  (2.14) 

n 

10'10 

180 

108 

— - —  — 1 

_q! 
a  10  *! 

- 

j 

V 

t 

j  2 

-oo  (l  +  ex) 

Eq.  (2.1) 

IT 

10"16 

168 

324 

— 

< 

¥ 


*A11  computations  were  carried  out  in  double  precision  floating  point  arithmetic 
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5.  PITFALLS  OF  CXXMPLTATION. 

The  accuracy  of  the  formulas  (2.1),  (2.6)  (2.10)  and  (2.14)  in  spite  of 
possible  singularities  at  the  end-points  of  an  interval,  is  based  on  our  being 
able  to  accurately  evaluate  the  integrand  in  a  neighborhood  of  these  singularities. 
This  is  especially  important  if  (as  is  often  the  case  for  singular  integrals) 
a  major  contribution  to  the  value  of  the  integral  occurs  in  a  neighborhood  of 
the  singularity. 

The  need  for  exercising  care  is  illustrated  by  considering  the  example 


(5.1) 


(3-2x-xV^dx 


:.hich  is  one  of  the  examples  in  Table  4.1.  Direct  application  of  Eq.  (2.14) 
to  the  approximation  of  I  results  in  the  formula 

L  kh  *  i 

(5.2)  Ia4h,  i  -fermf  (3-2^00-^00-1  * 

1~-K  (,e  + 1) 

where 


(5.3) 


.kh 


=  ~kh 


+  1 


The  points  z^Ch)  cluster  near  x  =  1  (resp.  near  x  =  -1)  for  k  large 

2 

and  positive  (resp.  large  and  negative).  Since  3-2x-x  =  0  when  x  *  1  , 
this  results  in  an  error  when  substituting  directly  into  (5.2)  to  evaluate 
this  quantity.  For  example,  if  h  =  \  ,  k  =  36  we  find,  working  to  8 
significant  figures,  that  z^(%)  ~  .9999  9997  ,  and  that 
[3-2zk(h)-z^(h)]”^  -•  (.0  v'0  001 2) =  2586. 7513.  On  the  other  hand,  if  x 


is 
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given  by  the  right  hand  side  of  (5.3),  we  get  [3-2x-x2}  ^  =  (%)  (e^1  +  l)(2e^t +1)  ^ 
=  2864.8728,  which  is  correct  to  8  significant  figures.  We  emphasize  that 
the  discrepancy  is  due  to  the  loss  of  significant  figures  in  the  evaluation  of 
z^(h)  via  the  use  of  (5.3). 

An  additional  error  occurs  in  the  evaluation  of  the  sun  in  (5.2).  In  terms 
of  z^(h)  ,  this  sum  may  be  written  in  the  form 

(5.2)'  I  a  2h  l  {l-z^(h)}[(3  +  zk(h))(l-zk(h))]'^  . 

k=-K 


The  trouble  occurs  in  the  computer  division  of  two  nuibers,  both  of  vfoich  are 
close  to  zero;  while  the  numerator  tern  {l-z£(h)}  expressed  in  the  form 
2ekh/[l  +  ekh]2  is  accurately  evaluated,  the  denominator  term  is  not,  since  the 
quantity  1  -  z,,(h)  only  has  1  significant  figure  of  accuracy. 

For  these  reasons  it  was  not  possible  to  achieve  an  error  <  10" 5  in  the 
evaluation  of  I  via  (5.2). 

If  we  replace  x  by  1  -  x  in  (5.1)  we  get  the  integral 


(5.4) 


dx 

[4x  -  x2]^ 


When  the  approximated  via  (2.14),  we  get  the  formula 


(5.5) 


J  2  4h 


L 

l 


to 


k=-K  (l  +  e^T 


[4xk(h)-x2(h)]"% 


where 


2e^ 

■  rl> 


(5.6) 
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The  integral  J  is  the  second  last  integral  in  Table  4.1.  The  singularity  of 
the  integrand  which  was  at  the  point  x  =  1  in  (5.1)  has  now  been  transformed 
to  the  point  x  =  0  in  (5.4).  In  contrast  to  the  loss  of  significant  figures 
encountered  in  the  evaluation  of  z^(h)  via  (5.3),  the  evaluation  of  z^Oi) 
via  (5.1)  can  be  carried  out  quite  accurately.  Nevertheless,  using  double 
precision,  we  were  still  only  able  to  achieve  10  significant  figures  of 
accuracy  in  the  approximation  of  J  via  (5.5).  The  reason  for  tl\is  is  the 
same  as  that  involving  the  discussion  of  (5.2)',  namely  requiring  the  computer 
to  evaluate  the  ratio  of  two  very  small  computed  quantities,  each  having  a 
slight  error. 

Finally,  let  us  replace  x  in  5.1  by  (ex-l)/(ex+l)  .  We  then  get  the 
integral 

,co  x 

(5.7)  H  S  2  — (l  +  2ex)"%dx  . 

•*-»  1  +  ex 

We  now  approximate  H  via  (2.1),  to  get 

i  -Wl  lyu  L 

(5.8)  H^2h  l  (1  +  26^)'^  . 

kr=-K  l  +  e™ 

In  the  expression  (5.7),  the  singularity  of  the  integrand  has  been  analytically 
removed.  We  thus  get  the  integral  in  the  last  entry  in  Table  4.1.  There  is 
now  no  problem  in  approximating  H  via  (5.8)  in  double  precision  to  get  16 
significant  figures  of  accuracy. 


- - 
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6.  SUBROUTINE  FOR  AUTOMATIC  INTEGRATION . 

The  FORTRAN  subroutine  for  evaluating  one  of  the  integrals  (1.1) -(1.4) 
is  called  by  the  statement 


WNEW(IOTRVL  ,EPS ,  A,  B ,  IP) 


Definitions  of  Parameters  and  Required  Function  Routine. 

(a)  General.  INTRVL,A,B  are  used  to  designate  the  domain  of  integration, 

and  which  integral  of  (1.1) -(1.4)  is  to  be  evaluated. 

EPS  is  a  small  positive  number,  specified  by  the  user.  This 
is  the  accuracy  to  which  the  integral  is  to  be  evaluated. 

IP  is  an  information  parameter. 

FN  is  the  name  of  a  user  supplied  function,  having  the 
form  FUNCTION  FN(X)  . 

(b)  Nbre  detailed  descriptions. 

INTRVL  =  1  means  that  the  integral  (1.1)  is  to  be  evaluated, 
user  sets  A  =  B  =  0  . 

INTRVL  =  2  means  that  the  integral  (1.2)  is  to  be  evaluated. 
User  sets  A  =  desired  numerical  value  as  in  (1.2),  B  «  0  . 

INTRVL  =  3  means  that  the  integral  (1 . 3)  is  to  be  evaluated. 
User  sets  A  =  desired  numerical  value  as  in  (1.3),  B  =  0  . 

INTRVL  =  4  means  that  the  integral  (1.4)  is  to  be  evaluated. 
User  defines  the  numerical  values  of  A  and  B  as  in  (1.4). 


IP  is  a  printout  information  parameter  selected  by  the  user.  It  is 
possible  to  have  the  following  lines  printed,  depending  on  values  of  IP 
(0,1  or  2)  chosen  by  the  user: 

(i)  H  LOWER  UPPER  T  M 

(ii)  D9.4  K  L  D20.10  D20.10 

(iii)  CONVERGENCE  ,  THE  FINAL  APPROXIMATION  IS  D30.17 

(iv)  DC  ENSIGN  EXCEEDED 

If  IP  =  0:  all  printouts  are  suppressed; 

IP  =  1 :  normal  printout  occurs .  This  includes  lines  (i)  (a  line  of  headings) , 
(ii)  (one  or  more  lines  of  numbers  under  the  headings  (i)),  and  one 
of  (iii)  or  (iv)  depending  on  whether  or  not  convergence  is 
achieved; 

IP  =  2:  only  line  (iv)  is  printed,  if  convergence  is  not  achieved. 

Let  us  briefly  explain  these  parameters  in  connection  with  what  the  program 
achieves . 

Let  us  denote  an  arbitrary  integral  (1.1) -(1.4)  by  I  . 

The  approximations  T  and  M  of  I  take  the  form 

L 

(6.1)  T  (=  Th)  -  h  l  w^(h)f(zk(h)) 

k=-K 

M  (=  ly  *  h  l  w2k-l^h/2)f(z2k-l(h^2))  • 


(6.2) 


The  nurerical  values  of  H  =  h  ,  K  ,  L  ,  T  and  M  in  line  (ii)  above  are 
the  parameters  in  (6.1)  and  (6.2).  The  integers  K  and  L  are  chosen  by 
the  program  (e.g.  for  T)  such  that 

-K+3 

(6.3)  l  Wjc(h)|f(zk(h))|  <  EPS/3  . 

k=-K 

L 

(6.4)  l  ^0016(^00)1  <  EPS/3  . 

k=L-  3 

Notice  that 

(6.5)  \n’^\+V  ■ 

Convergence  occurs,  and  the  printout  (iii)  follows  if  the  two  inequalities 

(6.6)  [T^l  <  EPS/3 
and 

(6.7)  K  +  L  +  1  <  5000 

are  satisfied.  In  this  case  the  nurber  given  by  (6.5)  is  printed  out 

in  line  (iii).  If  it  is  not  possible  to  achieve  the  requirements  (6.3),  (6.4) 
and  (6.6)  without  violating  (6.7),  the  error  message  (iv)  results. 

A  "sunrary"  flowchart  of  the  integration  routine  is  given  on  the  following 


page . 


! 


/  **  •  K 


SUBROUTINE  L'NEV  --  FLOWCHART 


H=1 

EPS3=F?S/3 

T 


CALL  GEM  —  Generate  nodes  x. (h)  and  corresponding  ; 

weights  w. (h)  for  the  upper  sum 

- -  -  i . -  - . 


Find  MAXUP  so  that  £j;j^,p_3wi(h)|f(xi(h))i  <  EPS3  : 
and  compute  U  “J.  w  (h)  f(x .(h)) 


w.(h)  f(Xi(h)) 


j  CALL  GEN  —  Generate  nodes  and  corresponding  weights  ! 

I  for  the  lower  sum  [ 

I*  ■  •  —  f 

!  Find  MAXLOW  so  that  y^AXL01',+w3.  (h)  (f(x.(h))|  <  EPS3 
I  ^i=MAXLO'r  1 

and  compute  L  *  Xi=MAXLOKwi  fO^Ch)) 

T  =  U  +  L  ' 

» 

I"  CALL  GENM  -  Generate  x^h/2)  and  corresponding  weights  w^(h/2) 

with  i  odd  and  satisfying  2  MAXLOW  +  1  5  i  <  2  MAXUP  -  1  : 


MAXUP-1 


M  =  h 


T  «21- 

i-MAXLOW 


3 (h/2)  f(x2il(h/2)) 


T  -  (T  +  >!)/2 


Yts 

\{T  -  Ml  EPS3  ■  ■*>- 


H  =  H/2 


>L\XUP  «  2  MAXUP 

MAXLOW  =  2  MAXI.OW 

* 

T  *  T 


PRINT  T 


STOP 

v/ 


'  ■* 
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7.  FORTRAN  LISTING  OF  WNEW. 

In  this  section  we  give  a  detailed  FORTRAN  listing  of  the  subroutine 
WNEW.  A  number  of  conments  are  given,  which  should  be  helpful  to  the  user 
of  this  routine. 


1 


nnoo  oooooooooooooo 


20 


SUBROUTINE  WNEW(INr,EPS,A,B,I?) 

C  UNEW  IS  THE  MAIN  SUBROUTINE  FOR  QUADRATURE . 


LNT-1  SIGNIFIES  AN  INFINITE  INTERVAL  -  THE  REAL  LINE . 

THEN  SET  A=B=0 

INT=2  AND  INT=3  SIGNIFIES  A  SEMI- INFINITE  INTERVAL 

(A.LT.X) .  THEN  SET  A  AS  THE  LEFT  HAND  END  POINT  AND 
B-0. 

INI-4  SIGNIFIES  A  FINITE  INTERVAL  (A.LT.X.LT.B) 

THEN  SET  A  AS  THE  LEFT  AND  B  AS  THE  RIQ{T  END 
POINT . 

EPS-DESIRED  DIFFERENCE  BETWEEN  THE  TRAPEZOIDAL  AND 
MIDPOINT  APPROXIMATIONS . 

FN-USER  SUPPLIED  FUNCTION  OF  THE  FORM  ’FACTION  F!(X)\ 
IP-0  SUPPRESSES  AIL  PRINTOUT. 

IP-1  FOR  NORMAL  PRE-TOUT. 

IP=2  FOR  ERROR  INDICATION  ONLY. 

WNEW  CALLS  SUBROUTINE  GEN  FOR  THE  INITIAL 
CALCULATION  OF  THE  TRAPEZOIDAL  APPROXIMATION.  AND 
THEREAFTER  CAII5  SUBROUTINE  GLNM  FOR  THE  APPROXIMATION 
OF  THE  1-TDPOLT  APPROXIMATION. 

IMPLICIT  DOUBLE  PRECIS  ION  (A-K.O-Z) 

DIMENSION  REIT (5000), VAL(5000) 

C  INITIALIZE  DATA. 

IF(IP.EQ.l)  TYPE  900 
H-1.0D0 
EPS3=EPS/3.0D0 
sun-o.  odo 

SUM2-0.0D0 

D=DSQRT(2.0D0) 

MAXP-48 

NBEG-1 

IF(INT.EQ. I) SUM=FN(0 . 0D0) 

IF(INT.EQ. 2) SUM=FN (A+I . 0D0) 

IF  (FT .  EQ .  3)  SLUHFN(A+DIOG  (L.ODOfD)  )  /D 
IF  (FT.  EQ .  4)  S0L?=  (B-A)  *FN  ( (A+B)  /2 .  ODO)  /4 . 0D0 


INITIALIZE  THE  UPPER  TAIL  OF  THE  TRAPEZOIDAL  APPROXIMATION 
CALL  GEN(INr,I,NBEC-.MA>T,A,B,WEIT,VAL) 

1=0 

CHEK=0.0D0 

mwxp 

DO  30  K=NBEG,M 
1=1+1 

EVA1M-EIT  (K)  *FN  ( VAL  (K) ) 

D=DABS(EVAL) 

CHEK=CHEK+D 

WEIT(K)=EVAL 

IF(I.LT.3)  GO  TO  30 

IF(CHEK.GT.EPS3)  GO  TO  20 

MAXP=K 

NBEG=1 

MAXL=MAX0 (MAXP , 48) 

GO  TO  35 
CHEK=0.0D0 
1=0 

CONTINUE 

NBEG=MAXP+1 

>:axp=maxp+48 

IF  (  (2,VMAXP)  .  GT .  5000)  GO  TO  110 
GO  TO  10 

DO  36  K=KAXP,1,-1 
SH-!2=Su-2H-EIT(K) 

INITIALIZE  THE  LOWER  TAIL  OF  THE  TRAPEZOIDAL 
APPROXIMATION. 

CALL  GEN(INr,2,NBEG,MAXL,A,B,WEIT,VAL) 

1=0 

CHEK=0.0D0 

M=MAXL 

DO  60  K=NBEG,M 
1=1+1 

EVALpWEIT  (K)  *FN  (VAL  (K)  ) 

WEIT(K)=EVAL 

D=DABS(EVAL) 

CHEK=CHEK+D 

IF(I.LT.3)  GO  TO  60 

IF(CHEK.GT.EPS3)  GO  TO  50 

MAXLfK 

GO  TO  70 

CHEKO.ODO 

1=0 

CONTINUE 

NBEG=MAXL+1 

MAXL^lAXL+48 

IF ( (MAXL-H 'AXP)  .  GT. 5000)  GO  TO  110 
GO  TO  40 

DO  75  K=MAXL,1,-1 
SUM1=SU!  T  fV7.  IT  (K) 

AP>T=Sl>l+5D!2+5UM 


COMPUTE  THE  MIDPOINT  APPROXTIMTION 
ML=MAXL 
12=>0.+1 

NTRM=MAXP4?AXL 

MAXLf=-MAXL 

NIRN*2*MAXL 

HOV2=K/2.0D0 

CALL  GE^T(IOT,OTRM,>rn?fLiOV2,H,A,B.WEIT,VAL) 

SUM=0 . 0D0 
SUML=O.ODO 
DO  85  K-l.Ml 

SUM=SLMWEIT(K)*FN(VAL(K)  ) 

DO  86  K=NTRM,M2,-1 
SUM=SUMI-H'.E  IT  (K)  *FN(VAL(K)  ) 
sum=silr5UM1 
AP»*-SIM*H 

TSTR-  (APXT+AP/P !)  /2 . 0D0 

IF(IP.EQ.  1)TY?E  901  ,H,MAXL,»:-?,A?}a ,APR»: 

IF  (DABS  (APXT-APXH)  .LT.EPS3)  C-0  TO  100 

SET  UP  DATA  FOR  THE  NEXT  ITERATION 

H=KOV2 

Xnw-2*NIRM 

M=MAXL 

M2«ML-I 

maxp^ax? 

APXT=TSTR 

IF(NTRM. GT.  5000)  GO  TO  110 
GO  TO  80 

REPORT  CONVERGENCE 
IF(IP.EQ.1)TYPE  902,  TSTR 
RETURN 

REPORT  FAILURE  TO  OBTAIN  CONVERGENCE 

IF(IP,GT.0)TYPE  903 

RETURN 

FORMAT (7X,  'H'  ,7X,  'LOWER'  ,3X,  ’UPPER'  ,11X,  ’T’  ,19X,  ’M’ ,/) 
FORMAT (3X , D9 , 4 , 218 , 2D20 . 10) 

FORMAT  (5X,  ’  CONVERGENCE  .THE  FINAL  APPROXIMATION  IS’,D30.I7) 
FOR?  AT ( 5X ,  ’ DIMENS ICWS  EXCEEDED') 

EJD 


SUBROUTINE  GEN(INT , LNF , NBEG , MAX , A , B .WEIT.VAL) 

IMPLICIT  DOUBLE  PRECISION (A-H.O-Z) 

DIMENSION  NEIT(5000) ,VAI,(5000) 

E=2 . 718281828459045235 36D0 

CALCULATION  of  THE  WEIGHTS  A'3  NODES  FOR  THE  TRAPEZOIDAL 
RULE 


INFINITE  INTERVAL 
IF (HJT.GE2)  GO  TO  10 
IF(INF.EQ.2)  GO  TO  6 
DO  5  K=N3EG,HAX 
WEIT(K)=1.0D0 
VAL  (K)  =DFL0AT  (K) 

CONTINUE 

RETURN 

DO  7  K=NBEG,MAX 

WEIT(K)=1.0D0 

VAL(K)=-DFLOAT(K) 

CONTINUE 
.  RETURN 

SEM  INFINITE  INTERVAL 
WEIT  (NBEG)  =E**N33EG 
DO  20  K=NBEC+1,K4X 
VEIT  (K)  =VJEIT  (K- 1)  *e 
CONTDXE 

IF(INF.EQ.l)  GO  TO  22 
DO  21  K=N3EG,!IAX 
VEIT(K)=1 .  QDO/WEIT  (K) 

IF(INT.EQ.3)GO  TO  26 

IF(INT.EQ.4)G0  TO  30 

DO  25  K=*3EG,tftX 

VAL(K)=A-..EITCQ 

CO'TTXE 

RETURN 

DO  29  K=NBEG,MAX 
W=WEIT(K) 

POt=DSQRT  (W)  *DSQRT  (1 . ODO/WW) 

IF(W.LT.O. 1D0)G0  TO  27 
VAL  (K)  =A+DL0G  C 0 
GO  TO  28 
W1=W*W 

W2- ( (-429 . DO/ 30720 .  DO*WI+231.  DO/13312.  DO)*Wl-63.  D0/2816.D0)*W1 

W2=  (  (  (W2+35 .  DO/ 1152 .  DO)  *Wl-5 .  DO/ 112 .  DO) *171+3 .  DO/40 .  DO)  *W1 

VAL(K)=(  (W2-1 .  DO/6 .  D0)*W1+1 .  DO)*W+A 

WEIT(K)=W/POM 

CONTINUE 

FINITE  INTERVAL 
5MA=B-A 

DO  40  K=NBEG,TAX 

DENM=WEIT(K)+1.0D0 

VAL  (K) = (  A+B'tyJE  IT  (K) )  /DENM 

weit(k)=b:'a*vjeit(k)  /  (de»we#o 

CONTINUE 

RETURN 
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SUBROUTINE  GENM(I\T , NTRM , NTRN ,  HOV2 ,  H , A ,  B , REIT , VAL) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DEENSION  VJEH\5000)  ,VAL(5000) 

C  CALCULATION  OF  THE  WEIGHTS  AND  NODES  FOR  THE  MIDPOINT  RULE 

C  LNFINITE  INVERAL 

IF(INT.GE.2)  C-0  TO  20 

NEIT(1)=1.0D0 

VAL(l)=DFLQAT(LHsTRN),,JHOV2 

DO  10  K=2,NTRM 

WEIT(K)=1.0D0 

VAL(K)-VAL(K-1)+H 

LO  CONTINUE 

RETURN 

C  SEMI  EJFINI7E  L'TERVAL 

20  EXFK=DEXP(ri) 

VEIT  (1)  =DEXP  (DFLOAI  (1+NT3N)  *H0V2) 

DO  30  K=2,NTRM 
VEIT(K)=WEIT  (K-IV'-EXPH 

30  CONTINUE 

IF(INT.EQ.  3)  GO  TO  36 
IF(L\T.EQ.4)  GO  TO  40 
DO  35  K=1,NTRM 
VAL  (K)  =A— ..7.  IT  (K) 

35  GQNTTNTE 
RETURN 

36  DO  39  K=1,NTPM 
W=WEIT(K) 

VJ1=W*W 

W3=DSQRT  (1 .  ODO-HvL) 

IFCW.LT.O. IDO)  GO  TO  37 
VAL(K)=DLOG(Vm%73)+A 
GO  TO  38 

37  W2=((-429 .  DO/  30720 .  DO*Wl+23l . 1X5/13312 .  D0),VW1~63 .  DO/2816 .  DO)  *W1 
W2=  (  (  (W2+35 .  DO/ 1 152 .  DO)  *Wl-5 .  DO/ 112 .  D0)*Wl+3 .  DO/40 .  D0)*W1 
VAL(K)=(  G'J2- 1 .  DO/  6 .  D0)*W1+1 . 0D0)*WA 

38  WEIT(K)=W/W3 

39  CONTINUE 
RETURN 

C  FINITE  INTERVAL 

40  BMA=B-A 

DO  50  K=lfNTRM 
DF.\M-WEIT(K)+1 . 0D0 
VAL(K)=(A+B"  (DM:.  -1 . 0D0) )  /DEM! 

WEIT  (K)  VV-  (DEM  I- 1 . 0D0)  /  (DKNM*DENM) 

50  CONTINUE 

RETURN 
ENT) 
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